# This script performs link-predictive analysis of the PE network
# Last edited: 9/24/2013
# Edited by: Bruce Desmarais

library(oc)

# Read in network data
load("~/Dropbox/professional/Research/Active/SenPE/Data/cleanData/netList.RData")
# amat is the adjacency matrix recording the number of events both attend
# vdat is the vertex dataset containing the variables in DSMKdata

### Cosponsorship ###
load("~/Dropbox/professional/Research/Active/SenPE/Data/cleanData/netListC.RData")

# functions to process vertex covariates
absdiffCov <- function(x){
	as.matrix(dist(cbind(x),diag=T,upper=T))
}

sameCov <- function(x){
	as.matrix(dist(cbind(x),diag=T,upper=T))==0
}


nodeCov <- function(x){
	matrix(x,length(x),length(x),byrow=T)+matrix(x,length(x),length(x),byrow=F)
}

overlapCov <- function(x,sep=":"){
	overMat <- matrix(0,length(x),length(x))
	for(i in 2:length(x)){
		for(j in 1:(i-1)){
			overMat[i,j] <- overMat[j,i] <- length(intersect(strsplit(x[i],split=sep)[[1]],strsplit(x[j],split=sep)[[1]]))
		}
	}
	overMat
}

listLengthCov <- function(x,sep=":"){
	xsp <- strsplit(x,split=sep)
	z <- unlist(lapply(xsp,length))
	nodeCov(z)
}

lowTri <- function(mat){
 	c(mat[lower.tri(mat)])
}

stDize <- function(mat){
	vecx <- c(mat)
	(mat - mean(vecx))/sd(vecx)
}



library(oc)
library(sna)
ests <- list()
set.seed(5)
for(i in 7:10){
	vdati <- netList[[i]]$vdat
	neti <- netList[[i]]$amat
	netCi <- netListC[[i]]$amat
	congi <- vdati$Congress[1]

	voteDat <- readKH(paste("ftp://voteview.com/sen",congi,"kh.ord",sep=""))

	icpsr <- voteDat$legis.data$icpsrLegis

	vdInNet <- match(vdati$ICPSR,icpsr)

	votes <- voteDat$votes
	votes[is.na(votes)] <- 0
	votes[votes==0] <- NA
	votes[is.element(votes,c(1,2,3))] <- 1
	votes[is.element(votes,c(7,8,9))] <- NA
	votes[is.element(votes,c(4,5,6))] <- -1
	
	votes <- votes[vdInNet,]

	corrYea <- cor(t(votes),use="pairwise")

	corrYea[corrYea==1] <- .99

	IdDist1 <- as.matrix(dist(vdati$DIm1,upper=T,diag=T))
	IdDist2 <- as.matrix(dist(vdati$Dim2,upper=T,diag=T))
	
	sameParty <- (as.matrix(dist(vdati$Party,upper=T,diag=T)) == 0)
	
	sameState <- (as.matrix(dist(vdati$StateCode,upper=T,diag=T)) == 0)

	X <- list(log(neti+1),log(netCi+1),IdDist1,IdDist2,sameParty,sameState)

	ests[[i]] <- netlm(qlogis((1+corrYea)/2),X,reps=1000)
	print(summary(ests[[i]]))
	print(i)
}

save(list="ests",file="~/Dropbox/professional/Research/Active/SenPE/Data/cleanData/RollCallResNL.RData")

# phi coefficient

congs <- c(96:105)

for(i in congs){
fileC <- paste("~/Dropbox/professional/Research/Active/SenPE/Tex/RCpredictC_",i,".pdf",sep="")
fileP <- paste("~/Dropbox/professional/Research/Active/SenPE/Tex/RCpredictP_",i,".pdf",sep="")

penet <- netList[[i-95]]$amat
cospnet <- netListC[[i-95]]$amat

xpe <- penet[lower.tri(penet)]
xcosp <- cospnet[lower.tri(cospnet)]
xsp <- seq(min(xpe),max(xpe),length=1000)
xsc <- seq(min(xcosp),max(xcosp),length=1000)

betap <- coef(ests[[i-95]])[2] 
betac <- coef(ests[[i-95]])[3]

xbp <- betap*log(xsp+1)
xbc <- betac*log(xsc+1)

sigp <- summary(ests[[i-95]])$pgreq[2] < .025
sigc <- summary(ests[[i-95]])$pgreq[3] < .025

yl <- c(min(c(xbp,xbc)),max(c(xbp,xbc)))

pdf(fileP,height=1.6, width=5.75, family="Times", pointsize=12.25)
par(las=1,mar=c(2.5,5.5,.75,.5),cex.lab=1.5,cex.axis=1.15)
plot(xsp,xbp,lty=ifelse(sigp,1,2),col=ifelse(sigp,"black","grey50"),lwd=4,type="l",xlab="",ylab="",ylim=yl)
#rug(xpe)
#title(xlab="Joint Press Events",line=2.25)
title(ylab="Corr",line=4)
dev.off()

pdf(fileC,height=1.6, width=5.75, family="Times", pointsize=11.25)
par(las=1,mar=c(2.5,5.5,.75,.5),cex.lab=1.5,cex.axis=1.15)
plot(xsc,xbc,lty=ifelse(sigc,1,2),col=ifelse(sigc,"black","grey50"),lwd=4,type="l",xlab="",ylab="",ylim=yl)
#rug(xcosp)
#title(xlab="Joint Cosponsorship",line=2.25)
title(ylab="Corr",line=4)
dev.off()

}






